import xarray as xr
from schimpy import schism_mesh
import pandas as pd
import holoviews as hv
from holoviews import opts, dim
ds = xr.open_mfdataset('../tests/data/m1_hello_schism/outputs/out2d_*.nc', concat_dim='time', combine="nested",
data_vars='minimal', coords='minimal', compat='override')
ds
<xarray.Dataset>
Dimensions: (time: 144, one: 1, nSCHISM_hgrid_node: 2639,
nSCHISM_hgrid_face: 4636,
nSCHISM_hgrid_edge: 7274,
nMaxSCHISM_hgrid_face_nodes: 4, two: 2)
Coordinates:
* time (time) datetime64[ns] 1999-12-31T16:20:00 ... 20...
SCHISM_hgrid_node_x (nSCHISM_hgrid_node) float64 dask.array<chunksize=(2639,), meta=np.ndarray>
SCHISM_hgrid_node_y (nSCHISM_hgrid_node) float64 dask.array<chunksize=(2639,), meta=np.ndarray>
SCHISM_hgrid_face_x (nSCHISM_hgrid_face) float64 dask.array<chunksize=(4636,), meta=np.ndarray>
SCHISM_hgrid_face_y (nSCHISM_hgrid_face) float64 dask.array<chunksize=(4636,), meta=np.ndarray>
SCHISM_hgrid_edge_x (nSCHISM_hgrid_edge) float64 dask.array<chunksize=(7274,), meta=np.ndarray>
SCHISM_hgrid_edge_y (nSCHISM_hgrid_edge) float64 dask.array<chunksize=(7274,), meta=np.ndarray>
Dimensions without coordinates: one, nSCHISM_hgrid_node, nSCHISM_hgrid_face,
nSCHISM_hgrid_edge,
nMaxSCHISM_hgrid_face_nodes, two
Data variables:
minimum_depth (one) float64 dask.array<chunksize=(1,), meta=np.ndarray>
SCHISM_hgrid (one) |S1 dask.array<chunksize=(1,), meta=np.ndarray>
depth (nSCHISM_hgrid_node) float32 dask.array<chunksize=(2639,), meta=np.ndarray>
bottom_index_node (nSCHISM_hgrid_node) int32 dask.array<chunksize=(2639,), meta=np.ndarray>
SCHISM_hgrid_face_nodes (nSCHISM_hgrid_face, nMaxSCHISM_hgrid_face_nodes) float64 dask.array<chunksize=(4636, 4), meta=np.ndarray>
SCHISM_hgrid_edge_nodes (nSCHISM_hgrid_edge, two) float64 dask.array<chunksize=(7274, 2), meta=np.ndarray>
dryFlagNode (time, nSCHISM_hgrid_node) float32 dask.array<chunksize=(72, 2639), meta=np.ndarray>
elevation (time, nSCHISM_hgrid_node) float32 dask.array<chunksize=(72, 2639), meta=np.ndarray>
depthAverageVelX (time, nSCHISM_hgrid_node) float32 dask.array<chunksize=(72, 2639), meta=np.ndarray>
depthAverageVelY (time, nSCHISM_hgrid_node) float32 dask.array<chunksize=(72, 2639), meta=np.ndarray>
dryFlagElement (time, nSCHISM_hgrid_face) float32 dask.array<chunksize=(72, 4636), meta=np.ndarray>
dryFlagSide (time, nSCHISM_hgrid_edge) float32 dask.array<chunksize=(72, 7274), meta=np.ndarray>This might not work for quad elements. Will have to deal with those as two non-overlapping triangles
smesh = schism_mesh.read_mesh('../tests/data/m1_hello_schism/hgrid.gr3')
dfelems = pd.DataFrame(smesh.elems,columns=[0,1,2])
#dfelems
dfnodes = pd.DataFrame(smesh.nodes, columns=['x','y','z'])
dfnodes.z = -dfnodes.z
#dfnodes
Needed to add simplices from existing elements information of mesh rather than Delaunay triangulation of nodes
The code below allows for the simplices already defined by elems to be used instead of doing a Delaunay triangulation (used from scipy as a way to calculate the simplices)
hv.extension('plotly')
import param
class TriSurface(hv.TriSurface):
simplices = param.Array()
class TriSurfacePlot(hv.plotting.plotly.TriSurfacePlot):
style_opts = ['cmap', 'plot_edges']
def get_data(self, element, ranges, style, **kwargs):
if element.simplices is None:
return super(TriSurfacePlot, self).get_data(element, ranges, style, **kwargs)
x, y, z = (element.dimension_values(i) for i in range(3))
simplices = element.simplices
return [dict(x=x, y=y, z=z, simplices=simplices)]
hv.Store.register({TriSurface: TriSurfacePlot}, 'plotly')
dfelev = dfnodes.copy()
#ds.elevation.min(), ds.elevation.max()
mesh_surface = TriSurface(dfelev, simplices = dfelems.values).opts(plot_edges=False, cmap='gray')
mesh_surface
The water surface simplices are derived from the nodes with Delaunay triangulation. Can we get the simplices from the information in the schism output files?
dfsurface = dfnodes.copy()
def show_surface(time=0):
dfsurface.z = ds.elevation.values[time,:]
water_surface = TriSurface(dfsurface, simplices = dfelems.values).opts(width=800, zlim=(-10,2))
return water_surface.opts(plot_edges=False, cmap='kbc', clim=(0,2), colorbar=True)
show_surface(0)
def show_combined(time):
return show_surface(time)*mesh_surface
import panel as pn
pn.extension()
time_slider = pn.widgets.IntSlider(name='Time Index', start=0, end=len(ds.time)-1)
pn.Row(time_slider,
hv.DynamicMap(show_combined, streams={'time':time_slider}).opts(
width=800, height=800)).servable(title='Water Surface Elevation Animation in 3D')#.show()